# Complementary material for submission to Transaction on Signal Processing
#
# Author: Wilder Lopes - wil@openga.org
# www.openga.org
# Feb 2018
#
# The following Python 3.0 (or above) script calls the GA-LMS (standard) binary to perform a system identification task.
#
# In the examples below, the multivectors samples belong to the Geometric Algebra of $\mathbb{R}^3$.
# Thus, each regressor and weight vector entry has 8 coefficients, i.e., for each entry of
# the weight vector, 8 coefficients have to be estimated. This constrasts with the usal LMS
# which only estimates real/complex entries. For further information, please refer to the
# GA documentation at www.openga.org.
#
# At the end, Fig. 5 (Top, Middle, and Bottom) and Fig. 6 of the submission
# are generated and saved in the folder "Figures".
#
#
# The MIT License (MIT)
#
# Copyright (c) 2016 onwards, by Wilder Lopes
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

# Start by importing the necessary Python modules:
import sys, string, os
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from math import log
from matplotlib.backends.backend_pdf import PdfPages

# Simulation parameters
#====================================================
# The user is able to set the following AF parameters:
#
# Number of filter taps (system order): M
#
# Realizations: L
#
# Time iterations: N
#
# AF step size: mu
#
# Measurement Noise variance: sigma2v

savefolder = './Figures/'
if not os.path.exists(savefolder):
    os.makedirs(savefolder)

###########################################
### Simulations to produce Fig. 5 (Top) ###
###########################################
print('### Experiments for Fig 5 (Top) ###')
pp = PdfPages(savefolder + '/TSP_GAAFs_Fig_5_Top.pdf') # multipage pdf to save figures

M       = 10 # System order
L       = 100 # Realizations
N       = 1000 # Time iterations
mu      = 0.01 # AF Step size
sigma2v = [1e-2, 1e-3, 1e-5] # Variance of measurement noise
sigma2q = 0 # Variance of random-walk noise
corr_input = 0.95 # Level of correlation between input's entries.
BINARY  = 'GA-LMS' # Note that you can call any of the following binaries:
		   # GA-LMS --> Complete subalgebra of R^3
		   # GA-LMS_rotors --> Even subalgebra of R^3 (isomorphic to quaternions)
		   # GA-LMS_complex --> Even subalgebra of R^2 (isomorphic to complex numbers)
		   # GA-LMS_real --> Even subalgebra of R (isomorphic to the real numbers)
#====================================================

data_dic = {}
for sig in sigma2v:
    # The binary is called below using the previously set parameters. The GA-LMS runs and returns .txt files with the results: *_galms.out and *_theory.out, where * represents "MSE" or "EMSE". *_galms.out files store the ensemble-average learning curves (EMSE), while *_theory.out files store the theoretical steady-state value for MSE and EMSE.

    # Calling binary
    arguments = " " + str(M) + " " + str(L) + " " + str(N) + " " + str(mu)+ " " + str(sig) + " " + str(sigma2q) + " " + str(corr_input)
    print('> Calling binary with the following parameters: \n {}'.format(arguments))
    os.popen("../../../src/GAAFs_standard/" + BINARY + "/build/" + BINARY + arguments).read()

    # Load files EMSE_galms.out and EMSE_theory.out (generated by the binary) to plot
    # EMSE learning curve and theoretical curve:
    f1 = open('EMSE_galms.out', 'r')
    data_label1 = ['EMSE_galms']
    data1_list = []
    for line in f1:
        data1_list.append(line.rstrip('\n'))

    f2 = open('EMSE_theory.out', 'r')
    data_label2 = ['EMSE_theory']
    data2_list = []
    for line in f2:
        for i in range(len(data1_list)):
            data2_list.append(line.rstrip('\n'))

    data1 = [float(j) for j in data1_list] # Converts to float
    data2 = [float(j) for j in data2_list] # Converts to float
    data1_dB = [10*log(x,10) for x in data1]
    data2_dB = [10*log(x,10) for x in data2]

    data_dic.update({'{}'.format(sig) : [data1_dB, data2_dB]})

# Plot figures
plt.title('EMSE curves - ' + BINARY + ' ' + r"$\mu={}, M={}$".format(mu, M))
plt.ylabel('EMSE (dB)')
plt.xlabel('Iterations')
for key in data_dic.keys():
    plt.plot(data_dic[key][1], label = r'$\sigma^2_v={}$ (theory)'.format(key), color = 'magenta', linestyle = '--')
    plt.plot(data_dic[key][0], label = r'$\sigma^2_v={}$'.format(key), color = 'b')

plt.legend(loc="best")
#plt.savefig('EMSE.png', bbox_inches='tight')
pp.savefig()
#plt.show()
plt.close()
pp.close()



##############################################
### Simulations to produce Fig. 5 (Middle) ###
##############################################
print('### Experiments for Fig 5 (Middle) ###')
pp = PdfPages(savefolder + 'TSP_GAAFs_Fig_5_Middle.pdf') # multipage pdf to save figures

M       = 10 # System order
L       = 100 # Realizations
N       = 1000 # Time iterations
mu      = [0.005, 0.0075, 0.01, 0.0125, 0.015, 0.0175, 0.02, 0.0225] # AF Step size
sigma2v = [1e-2, 1e-3, 1e-5] # Variance of measurement noise
sigma2q = 0 # Variance of random-walk noise
corr_input = 0.95 # Level of correlation between input's entries.
BINARY  = 'GA-LMS' # Note that you can call any of the following binaries:
		   # GA-LMS --> Complete subalgebra of R^3
		   # GA-LMS_rotors --> Even subalgebra of R^3 (isomorphic to quaternions)
		   # GA-LMS_complex --> Even subalgebra of R^2 (isomorphic to complex numbers)
		   # GA-LMS_real --> Even subalgebra of R (isomorphic to the real numbers)
#====================================================

# Number of points to average the EMSE:
avgpoints = 200

data_dic = {}
for sig in sigma2v:
    data1_dB = []
    data2_dB = []
    for step in mu:

        # Calling binary
        arguments = " " + str(M) + " " + str(L) + " " + str(N) + " " + str(step)+ " " + str(sig)  + " " + str(sigma2q) + " " + str(corr_input)
        print('> Calling binary with the following parameters: \n {}'.format(arguments))
        os.popen("../../../src/GAAFs_standard/" + BINARY + "/build/" + BINARY + arguments).read()

        # Load files EMSE_galms.out and EMSE_theory.out (generated by the binary) to plot
        # EMSE learning curve and theoretical curve:
        f1 = open('EMSE_galms.out', 'r')
        data_label1 = ['EMSE_galms']
        data1_list = []
        for line in f1:
            data1_list.append(line.rstrip('\n'))

        f2 = open('EMSE_theory.out', 'r')
        data_label2 = ['EMSE_theory']
        data2_list = []
        for line in f2:
            for i in range(len(data1_list)):
                data2_list.append(line.rstrip('\n'))

        data1 = [float(j) for j in data1_list] # Converts to float
        data2 = [float(j) for j in data2_list] # Converts to float

        data1_dB.append(10*log(np.mean(data1[-avgpoints:]),10))
        try:
            data2_dB.append(10*log(data2[0],10))
        except:
            data2_dB.append(np.nan)

    data_dic.update({'{}'.format(sig) : [data1_dB, data2_dB]})

# Plot figures
plt.title('EMSE curves - ' + BINARY + ' ' + r"$M={}$".format(M))
plt.ylabel('EMSE (dB)')
plt.xlabel('Step size')
for key in data_dic.keys():
    plt.plot(mu, data_dic[key][1], label = r'$\sigma^2_v={}$ (theory)'.format(key), color = 'magenta', linestyle = '--', marker='o')
    plt.plot(mu, data_dic[key][0], label = r'$\sigma^2_v={}$'.format(key), color = 'black', marker='x')

plt.legend(loc="best")
#plt.savefig('EMSE.png', bbox_inches='tight')
pp.savefig()
#plt.show()
plt.close()
pp.close()


##############################################
### Simulations to produce Fig. 5 (Bottom) ###
##############################################
print('### Experiments for Fig 5 (Bottom) ###')
pp = PdfPages(savefolder + 'TSP_GAAFs_Fig_5_Bottom.pdf') # multipage pdf to save figures

M       = [x for x in range(1,25)] # System order
L       = 100 # Realizations
N       = 1000 # Time iterations
mu      = 0.01 # AF Step size
sigma2v = [1e-2, 1e-3, 1e-5] # Variance of measurement noise
sigma2q = 0 # Variance of random-walk noise
corr_input = 0.95 # Level of correlation between input's entries.
BINARY  = 'GA-LMS' # Note that you can call any of the following binaries:
		   # GA-LMS --> Complete subalgebra of R^3
		   # GA-LMS_rotors --> Even subalgebra of R^3 (isomorphic to quaternions)
		   # GA-LMS_complex --> Even subalgebra of R^2 (isomorphic to complex numbers)
		   # GA-LMS_real --> Even subalgebra of R (isomorphic to the real numbers)
#====================================================

# Number of points to average the EMSE:
avgpoints = 200

data_dic = {}
for sig in sigma2v:
    data1_dB = []
    data2_dB = []
    for order in M:

        # Calling binary
        arguments = " " + str(order) + " " + str(L) + " " + str(N) + " " + str(mu)+ " " + str(sig)  + " " + str(sigma2q) + " " + str(corr_input)
        print('> Calling binary with the following parameters: \n {}'.format(arguments))
        os.popen("../../../src/GAAFs_standard/" + BINARY + "/build/" + BINARY + arguments).read()

        # Load files EMSE_galms.out and EMSE_theory.out (generated by the binary) to plot
        # EMSE learning curve and theoretical curve:
        f1 = open('EMSE_galms.out', 'r')
        data_label1 = ['EMSE_galms']
        data1_list = []
        for line in f1:
            data1_list.append(line.rstrip('\n'))

        f2 = open('EMSE_theory.out', 'r')
        data_label2 = ['EMSE_theory']
        data2_list = []
        for line in f2:
            for i in range(len(data1_list)):
                data2_list.append(line.rstrip('\n'))

        data1 = [float(j) for j in data1_list] # Converts to float
        data2 = [float(j) for j in data2_list] # Converts to float

        data1_dB.append(10*log(np.mean(data1[-avgpoints:]),10))
        try:
            data2_dB.append(10*log(data2[0],10))
        except:
            data2_dB.append(np.nan)

    data_dic.update({'{}'.format(sig) : [data1_dB, data2_dB]})

# Plot figures
plt.title('EMSE curves - ' + BINARY + ' ' + r"$\mu={}$".format(mu))
plt.ylabel('EMSE (dB)')
plt.xlabel('System Order M (Taps)')
for key in data_dic.keys():
    plt.plot(M, data_dic[key][1], label = r'$\sigma^2_v={}$ (theory)'.format(key), color = 'magenta', linestyle = '--', marker='o')
    plt.plot(M, data_dic[key][0], label = r'$\sigma^2_v={}$'.format(key), color = 'black', marker='x')

plt.legend(loc="best")
#plt.savefig('EMSE.png', bbox_inches='tight')
pp.savefig()
#plt.show()
plt.close()
pp.close()


##############################################
###      Simulations to produce Fig. 6     ###
##############################################
print('### Experiments for Fig 6 ###')
pp = PdfPages(savefolder + 'TSP_GAAFs_Fig_6.pdf') # multipage pdf to save figures

M       = 10 # System order
L       = 100 # Realizations
N       = 2000 # Time iterations
mu      = 0.005 # AF Step size
sigma2v = 1e-3 # Variance of measurement noise
sigma2q = 0 # Variance of random-walk noise
corr_input = 0.95 # Level of correlation between input's entries.
BINARY  = ['GA-LMS_rotors', 'GA-LMS_complex', 'GA-LMS_real']
		   # GA-LMS --> Complete subalgebra of R^3
		   # GA-LMS_rotors --> Even subalgebra of R^3 (isomorphic to quaternions)
		   # GA-LMS_complex --> Even subalgebra of R^2 (isomorphic to complex numbers)
		   # GA-LMS_real --> Even subalgebra of R (isomorphic to the real numbers)
#====================================================

# Number of points to average the EMSE:
avgpoints = 200

color_dic = {'GA-LMS_rotors' : 'blue',
             'GA-LMS_complex' : 'green',
             'GA-LMS_real' : 'black'}
data_dic = {}
for binary in BINARY:
    data1_dB = []
    data2_dB = []

    # Calling binary
    arguments = " " + str(M) + " " + str(L) + " " + str(N) + " " + str(mu)+ " " + str(sigma2v)  + " " + str(sigma2q) + " " + str(corr_input)
    print('> Calling binary with the following parameters: \n {}'.format(arguments))
    os.popen("../../../src/GAAFs_standard/" + binary + "/build/" + binary + arguments).read()

    # Load files EMSE_galms.out and EMSE_theory.out (generated by the binary) to plot
    # EMSE learning curve and theoretical curve:
    f1 = open('EMSE_galms.out', 'r')
    data_label1 = ['EMSE_galms']
    data1_list = []
    for line in f1:
        data1_list.append(line.rstrip('\n'))

    f2 = open('EMSE_theory.out', 'r')
    data_label2 = ['EMSE_theory']
    data2_list = []
    for line in f2:
        for i in range(len(data1_list)):
            data2_list.append(line.rstrip('\n'))

    data1 = [float(j) for j in data1_list] # Converts to float
    data2 = [float(j) for j in data2_list] # Converts to float
    data1_dB = [10*log(x,10) for x in data1]
    data2_dB = [10*log(x,10) for x in data2]

    data_dic.update({'{}'.format(binary) : [data1_dB, data2_dB]})


# Plot figures
plt.title('EMSE curves - ' + r"$M={}, \mu={}, \sigma^2_v={}$".format(M, mu, sigma2v))
plt.ylabel('EMSE (dB)')
plt.xlabel('Iterations')
for key in data_dic.keys():
    plt.plot(data_dic[key][0], label = '{}'.format(key), color = color_dic[key])
    plt.plot(data_dic[key][1], label = '{} (theory)'.format(key), color = 'magenta', linestyle = '--')

plt.legend(loc="best")
#plt.savefig('EMSE.png', bbox_inches='tight')
pp.savefig()
#plt.show()
plt.close()
pp.close()
